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ERROR ESTIMATION AND UNCERTAINTY PROPAGATION 
IN COMPUTATIONAL FLUID MECHANICS 

J.Z. ZHU* AND GUOWEI H F. t 

Abstract. Numerical simulation has now become an integral part of engineering design process. Critical 
design decisions are routinely made based on the simulation results and conclusions. Verification and vali- 
dation of the reliability of the numerical simulation is therefore vitally important in the engineering design 
processes. We propose to develop theories and methodologies that can automatically provide quantitative 
information about the reliability of the numerical simulation by estimating numerical approximation error, 
computational model induced errors and the uncertainties contained in the mathematical models so that 
the reliability of the numerical simulation can be verified and validated. We also propose to develop and 
implement methodologies and techniques that can control the error and uncertainty during the numerical 
simulation so that the reliability of the numerical simulation can be improved. 

Key words, error and uncertainty in numerical simulation, ZZ error estimators, error energy spectrum, 
mapping closure approximation, adaptive mesh refinement 

Subject classification. Fluid Mechanics 

1. Introduction. Our ultimate goal of the research is to develop a computation-based procedure that is 
capable of automatically estimating both numerical and model-induced errors as well as controlling the level 
of uncertainty in physical models, so that the physical quantities of engineering interests can be correctly 
predicted. The research is therefore falling into the scope of verification and validation of the reliability 
of the numerical simulation, with emphases on both computational solid mechanics and computational 
fluid dynamics (CFD). There are two distinct issues associated with verification and validation of numerical 
simulation, i.e., control of the error in numerical computations and computational models, and quantification 
of the uncertainty in physical models. Here error is considered as a recognizable deficiency in any phase of 
simulation that is not due to lack of knowledge. Uncertainty is, on the other hand, defined as a potential 
deficiency in any phase or activity of the modeling and simulation that is due to lack of knowledge. 

Numerical simulation has now become an integral part of engineering design process. Critical design 
decisions are routinely made based on the simulation results and conclusions. An inaccurate design due 
to unreliable numerical simulation, caused by error and uncertainty, can result in waste, performance loss 
and even catastrophe. This is particularly true in aerospace and automotive industries. To achieve reliable 
numerical simulations that can be confidently used in engineering design process, quantitative measurement 
and effective control of the error and uncertainty are vitally important. By accurately assessing the error 
appeared in the numerical computations, computational models, and uncertainties contained in the physical 
models, we will be able to develop numerical algorithms that not only improve the accuracy of the numerical 
simulation, but also provide users with the knowledge of the reliability of the numerical results. 

In what follows, we shall introduce the state of the art, which includes our own efforts, of a posteriori 
error estimation of the numerical simulation, error analysis of the computational models and uncertainty 
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propagation in physical models. A posteriori error estimation of the numerical approximation, the essential 
ingredient of the verification of the numerical simulation, is addressed in Section 2. The possibility of error 
analysis in computational models is discussed in Section 3. A methodology is proposed, in Section 4, to 
quantify the uncertainty propagation in physical models. Section 5 describes an adaptive computational 
method, which is an optimal tool to control the accuracy of the numerical solution. Finally, in Section 6, 
the proposed research and development are presented. 

2. A posteriori error estimation. Mastering the numerical modeling to achieve an accurate and 
reliable solution for a particular engineering problem described by certain mathematical models is now a 
crucial part of the modern computational mechanics. The numerical method, such as the finite element 
method or finite difference method, leads to approximate solutions of the mathematical model; the quality 
of these solutions must be verified so that the numerical predictions can be validated. A posteriori error 
estimation plays a central role in the verification of the reliability of the approximate solutions. 

Error estimations are called as a posteriori when they are constructed after the approximate solution has 
been computed. A posteriori error estimation utilize the numerical solution and sometimes the input data so 
that the actual error presented in the approximate solution can be calculated. In numerical computations, 
the numerical approximation error is a function of position and time and is defined as eh =u-Uh, where u 
is the exact solution of the mathematical model that represents the physical phenomenon of interest; Uh is 
the approximate solution of certain numerical method, such as a finite element method. The error can be 
measured and estimated under appropriate norms. The estimated error is called a posteriori error estimator. 

There are essentially two types of error estimators: 

(1) The residual type of error estimators. This type of error estimators, first introduced by Babuska 
and Rheinboldt [1] and later further developed by Zienkiewicz et al. [2], Bank and Weiser [3], Ainsworth 
and Oden [4], and others [5], is computed by using the residuals of the approximate solution, which fail to 
satisfy the mathematical model, explicitly or implicitly. The explicit residual error estimators are expressed 
directly by the residuals of the numerical solution; the implicit residual error estimators, on the other hand, 
are determined by solving local equations. The explicit residual error estimators have been show to be not 
accurate in estimating the error, and can only be used for guiding the adaptive analysis procedures [6, 7]. 
The implicit residual error estimators are more accurate, but needs substantial computational effort; the 
requirement of solving local problem also makes it difficult to be applied to non-linear problems. 

(2) The recovery type of error estimators. This type of error estimators, proposed by Zienkiewicz and 
Zhu [8] and often referred as ZZ error estimators, is now widely used by researchers and adopted by almost 
all the major commercial numerical modeling software. The error estimators are calculated by first using 
some recovery techniques to achieve a more accurate, or even superconvergent solution from the computed 
numerical solutions. The recovered solution is then used in place of the exact solution to compute the 
error. When the superconvergent recovery (SPR) technique [9, 10, 11] is used as the recovery technique, this 
type of error estimators is remarkably accurate and extremely robust [6, 7]. The methodology of this type 
of error estimators is quite general, although it was first intended for finite element methods. A recovery 
error estimator is independent of the mathematical models being solved, therefore applicable to non-linear 
problems [12, 13]; it is also independent of the discretization method, therefore applicable not only to finite 
element methods but also finite difference and finite volume methods [14]. This type of error estimators is 
computationally more efficient. 

The above error estimators, in their current forms, only give estimation of the error in certain norms; it 
does not directly provide information about the error for a particular quantity of interest, for example, the 
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error in temperature, heat flux, velocity, pressure, displacements or stresses etc. To effectively compute the 
point-wise error of the quantity of interest, a recovery type of point-wise error estimator has to be developed. 

The vital role that a posteriori error estimation plays in assessing the accuracy of the numerical solution 
approximating to the mathematical model is well recognized. It is at the center of the verification of the 
reliability of the numerical simulation. A posteriori error estimator is also essential for adaptive analysis 
procedure which controls the error within a prescribed tolerance by systematically and automatically refining 
the mesh, the numerical algorithms and even the discretization method. 

3. Error estimation in computational models. One of the goals of our research is to develop a 
computation-based procedure which can automatically select appropriate turbulence models used in numer- 
ical simulation that can correctly predict the physical quantities of engineering interests. This will be based 
on the estimation of the errors induced by turbulence models. In order to estimate the error of the turbulence 
model properly, it is assumed that the discretization errors are not dominant in the computational results 
so that the effects of the model errors can be isolated. The discretization errors have been discussed in the 
preceding section. 

Turbulence models are mathematical models for turbulence flow. The mathematical models used in 
numerical simulations are usually referred to as computational models. Therefore, we will address the issue 
of validation of the reliability of the computational models used in the numerical simulation and develop 
ways to hierarchically improve the mathematical models when the accuracy of these models is found not to 
be sufficient. The computational models of a physical problem are derived under various assumptions and 
simplifications. For the same physical problem, different levels of assumption and simplification result in a 
hierarchy of computational models. This can been seen in a variety of turbulence models similar to that 
appeared in solid mechanics in forms of plate bending formulation. The adoption of a much simplified com- 
putational model in solving the physical phenomenon of interest is often based on heuristics and experience 
in hope that the problem can be solved effectively with sufficient accuracy. The numerical solutions of such 
simplified computational models are often neither validated with the more sophisticated models, nor verified 
with the physical phenomena. The conclusions derived from the comparison of the numerical solutions to 
these models with the equally simplified experiments can sometimes be misleading. The estimation of the 
errors induced by the simplified computational model is, therefore, particularly important in determining 
the accuracy of the numerical simulation for the practical problem considered. 

The theory and technique of the error estimation on the computational models are much less developed in 
contrast to the accuracy assessment of the numerical solution for a given mathematical model [15] [16]. This 
is mainly due to the fact that the exact solutions of physical problems are usually unknown and the direct 
numerical simulation (DNS) solutions are practically inaccessible because of the enormous computational 
resource requirement. Thus, accurate and efficient a posteriori error estimation of the computational models 
is highly desirable in CFD. Some preliminary research work has been done in the adaptive hierarchic modeling 
for computational models, in the context of solid mechanics, by Babuska et al [17], Oden et al [18] and 
Ainsworth et al [19] in the last ten years; Much research and development are needed in CFD and particularly 
in the turbulence flow computation. In the following, we wall explore possible theories and techniques of 
error estimation on turbulence modeling. 

(1) Direct error estimation 

The error estimation on turbulence modeling can be directly achieved by comparing the results from a 
simulation that uses the turbulence models with the ones from DNS, if the DNS data, such as mean velocity 
and Reynolds stress distributions, are available. Their differences provide an estimation of the error induced 
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by the turbulence models. The estimation can be considered as exact if the DNS is considered as exact. 
Unfortunately, the DNS data are usually not available for practical problems. The direct error estimation 
is. therefore, practically impossible in CFD. 

(2) Indirect error estimation 

The basic idea of the indirect error estimation is to find a surrogate of the DNS solution. The difference 
between numerical solution with turbulence model and the surrogate provides an error estimation of the 
turbulence model. The surrogate may be found by either solving surrogate equations or using extrapolation 
techniques. For example, a more sophisticated computational model, instead of the exact mathematical 
model such as Navier-Stokes equations, could be solved for a quality assessment of the simplified models, 
although this will require overwhelming computational time and cost. 

We are working on development of extrapolation techniques that can utilize multi-level turbulence models 
to quantitatively assess the accuracy of the turbulence models. The possibility of using a technique analogous 
to the one used in the estimation of the discretization error for a given mathematical model, as these 
discussed in Section 2, will be exploited. The main idea of developing a feasible method is to use solutions 
for turbulence models at different levels most economically together with a proper extrapolation technique 
to estimate the error of a particular computational model and its numerical approximation. We will develop 
the extrapolation techniques that can use numerical solutions from different mathematical models to produce 
a solution with the same accuracy as the most complicated model. 

(3) The statistical estimates: error energy spectrum 

Although it is impossible in deterministic fashion, the direct estimation of the error may be achieved 
in statistical fashion, such as moments and Probability Density Functions (PDF). The moment and PDF 
equations can be derived and their solutions provide the probabilistic measurement of the error. For example, 
the second-order moments of the error are the standard deviation of the model solutions from the exact 
solutions. The energy spectra of the error describe the distribution of total error at different scales. The 
second author’s recent research [20] indicates that the model error could be amplified but bounded by the 
total energy. The error amplification is accompanied by an inverse cascade process of errors from small scales 
to large scales. 

The concept of the error energy spectrum was initially developed for the predictability problem of 
turbulence, i.e., the propagation of the initial errors in Navier-Stokes equation [21, 22], Here we suggest to 
use it in the analysis of the error induced by turbulence models. The mathematical equation of the error 
energy spectrum can be formulated from the Navier-Stokes equation and the model equations, and it is 
closed by analytical theory of turbulence. The equation of error energy spectrum is much easier to solve 
than those equations involved in the direct and indirect error estimation. Analytical approximations that 
include statistics no higher than two-time, two-point velocity variance has afforded adequate prediction of 
the evolution of turbulence spectra [23]. Similar methods have been succeeded in predicting the growth of 
initial errors. Therefore, it is expected that the error energy spectrum approach will provide a satisfactory 
estimation of the model-induced errors. 

4. The propagation of uncertainty in physical models. Uncertainties are inevitably contained 
in physical models due to lack of knowledge [15, 16]. These uncertainties might be present in boundary 
conditions, initial conditions and model parameters. How various uncertainties affect the reliability of the 
numerical simulation is of major concern in engineering design. The influence of the uncertainty in the 
physical model on the outcome of the numerical simulation is referred as uncertainty propagation. 

Uncertainty propagation appeals for a posteriori probabilistic description of the numerical simulation. 


4 



This usually requires either moment or PDF approach. However, both moment and PDF approaches suffer 
from closure problems, i.e., there are more unknown terms than the available transport equations. Those 
additional unknowns have to be approximated by certain models so that the problem can be solved, as 
appeared in turbulence modeling where the closure problems are often addressed by using Kolmogorov 
theory. Unfortunately, such a sound theory does not exist in quantifying the uncertainty propagation. 

An effective approach to uncertainty propagation is the nonlinear mapping approximation [24]. An 
augment mapping [25] from initial conditions and/or random input to random output is defined through 
governing equations. The statistics of the output can be calculated as soon as the mapping is known. The 
transport equations of the nonlinear mapping are formulated from the governing equations and its PDF 
equations. 

A direct application of the augment mapping is to solve the governing equation for all realizations of 
the input and compute all realizations of the output. An average of all realizations of the output gives the 
desired statistics. This approach is similar to the Monte Carlo method, but it is prohibitively expensive. 

The Taylor expansion of the mapping with second order truncation is probably the easiest way to 
calculate the first and second order moments of the output [26]. But it works only for small variance due to 
approximation errors. 

The transport equation of the nonlinear mapping does not, in fact, need to be solved exactly, since the 
mapping is only required to have the same statistics as the exact solution of the governing equation. In other 
words, we require that the PDF of the nonlinear mapping obey the PDF transport equation of the output. 
Based on the PDF equation, it is feasible to formulate a governing equation for the mapping. 

As the mapping equation is in the closed form that does not need any ad hoc models, it constitutes 
a primary advantage for the mapping approximation over moment and PDF approaches where modeling 
unclosed terms is unavoidable. The mapping equation is deterministic because it also avoids Monte Carlo 
simulation. The mapping approximation can be performed numerically by a numerical method such as finite 
difference method. An application of the augment mapping approach to quasi 1-D Euler CFD problem shows 
that it gives better prediction of uncertainty propagation [27]. 

5. Adaptive computational methods. Adaptive computational methods refer to algorithms for var- 
ious numerical methods, for example, finite element methods, finite difference methods and finite volume 
methods, which allow the size and the distribution of the computational mesh, the order of the approxi- 
mation and the time steps to be changed dynamically, guided by a posteriori error estimators, so that the 
discretization error can be controlled within the given tolerance. In the context of the finite element method, 
the modification of the mesh size and mesh distribution is termed as h- version adaptive refinement, where 
the error is reduced by adaptively refining the mesh and leaving the order of the approximation unchanged; 
On the other hand, the modification of the approximation order is named as p- version adaptive refinement, 
where the error is reduced by increasing the order of the approximation function while keeping the mesh 
fixed. A proper combination of h-version and p-version is the hp-version adaptive refinement. The algo- 
rithms that involve mesh refinement are more challenging than a mere order increasing. The spatial mesh 
refinement, however, poses much more challenges than the modification of time steps. These are particularly 
true for 3-D problems. Therefore, we shall concentrate mainly our research and development in the spatial 
mesh refinement, i.e., h-version mesh refinement, and its related problems that appear in 3-D problems and 
in non-linear problems. 

There are two types of mesh refinement techniques currently used in h-version adaptive computations: 
(1) Local mesh refinement (or mesh enrichment); (2) Re-meshing (or mesh regeneration). With local mesh 
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refinement the new mesh is obtained by bisecting the mesh locally so that the elements with the largest 
local estimated error are refined, based on some appropriate refinement criteria at each step of the adaptive 
analysis procedure. The main feature of the local mesh refinement is that it allows the error distribution to 
be controlled precisely. In the process of the mesh refinement, only the elements where the errors exceed the 
prescribed level will be refined. The resulting mesh of this process is in general very close to the optimal 
mesh, i.e., a mesh with minimum number of elements for a specified accuracy. 

The advantages of using local mesh refinement are that, for a given geometry and an initial mesh, the 
mesh refinement can be carried out rather easily as long as the data structure that manages the relationship 
of the elements from different levels of refinement is in place. However, for problems that involve a change 
of the geometry, such as simulations of hypervelocity impact, optimum design and metal forming, local or 
often complete mesh regeneration is inevitable. Adaptive re-meshing will be the natural choice to be used 
in numerical solution of these kinds of problems. 

Adaptive re-meshing depends on not only the reliability of a posteriori error estimator, but also the 
availability of an automatic mesh generator that can, following a proper mesh refinement strategy, produce 
meshes with nearly optimal mesh size distribution for a prescribed accuracy. 

Although both local mesh refinement and re-meshing have long been used in the numerical simulations 
benefited from the appearance of the automatic mesh generators [28]-[33], some issues related with non-linear 
adaptive computation still need to be resolved properly. One of the main issues is the proper transfer (or 
mapping) of the data, such as history dependent parameters and solid-fluid interfaces etc., when elements 
are refined. The robustness of adaptive analysis for non-linear problems is in fact largely depends on the 
correctness of the information transferred from an old mesh to a new mesh. Most of the methodologies 
suggested in the literature are using certain type of interpolation techniques in transferring data [34]-[36], 
however, none of them seems to be totally robust. 

The availability of the automatic mesh generators that can generate meshes predicted by error estimator 
in both 2-D and 3-D is another issue. Currently, few 3-D tetrahedral mesh generators have such capability. 
Automatic hexahedral mesh generator does not even exist. Both tetrahedral mesh generation and hexahe- 
dral mesh generation have been active research subjects heavily funded by the United States government, 
European Union and private companies in the last 20 years. It is one of the most competitive research areas 
in numerical simulations both academically and commercially, because of its tremendous potential academic 
and commercial values. 

6. Future work. The reliability of the numerical simulation is critical in the decision-making of engi- 
neering design. The error that created by numerical approximation, induced by mathematical and physical 
models, if not properly estimated and effectively controlled, can render the solutions of the numerical sim- 
ulation totally invalid. Design decisions based on such solutions, which almost always lead to waste, can 
sometimes lead to disastrous consequences. The verification of the numerical computation and validation of 
the mathematical and physical models are vitally important in any numerical simulation. The objective of 
this research is to develop an integrated approach to verify the numerical simulation by a posteriori error 
estimation of the numerical solutions, validate the mathematical model and physical model by developing 
efficient technique including nonlinear mapping approximation to estimate the model induced error, and 
effectively reduce errors by using adaptive computational method for engineering numerical simulations. 
Utilizing local mesh refinement technique, together with automatic selection of efficient computational mod- 
els and control of uncertainty in physical models will be able to create credible numerical simulations for 
engineering design. 
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The expected results that will be produced as the result of this research include: 

(1) The recovery type of the error estimators will be implemented and its extension for the point-wise 
error estimation of the quantity of interest will be developed for problems of engineering applications, to 
provide verification for the reliability of the numerical simulations. The recovery type of the error estimators 
will be incorporated in the adaptive computational method. 

(2) Extrapolation techniques that can use numerical solutions from different mathematical models for 
error estimation of the computational models will be developed. Statistical approaches including error energy 
spectrum method will be also developed. 

(3) The nonlinear-mapping-approximation based probabilistic approach for uncertainty propagation will 
be developed, which can be used to calculate moment and probability density function of uncertainty in 
numerical simulations. 

(4) The existing techniques for automatic triangular and quadrilateral mesh generators, which are capable 
of generating meshes suitable for adaptive computational methods, in 2-D and tetrahedral mesh generator 
in 3-D will be implemented. Research and development of hexahedral mesh generator will be conducted. 

(5) New techniques for robust information transfer will be developed. Here, superconvergence recovery 
technique can be applied. 
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